Simulation of Unsaturated Flow in Complex Fractures Using Smoothed Particle Hydrodynamics

on SPH to simulate three-dimensional invasion of nonwetting and wetting fluids in fracture apertures with Smoothed particle hydrodynamics (SPH) models were used to complex geometries. Such complex processes are very simulate unsaturated flow in fractures with complex geometries. The SPH is a fully Lagrangian particle-based method that allows the dydifficult to model with traditional grid-based numerical namics of interfaces separating fluids to be modeled without emschemes. The SPH is an interpolation-based technique ploying complex front tracking schemes. In SPH simulations, the fluid that can be used to solve numerically systems of partial density field is represented by a superposition of weighting functions differential equations. The Lagrangian particle nature centered on particles which represent the fluids. The pressure is related of SPH allows physical and chemical effects to be incorto the fluid density through an equation of state, and the particles porated into the modeling of flow processes, which may move in response to the pressure gradient. The SPH does not require also be complicated by irregular and deformable boundthe construction of grids that would otherwise introduce numerical aries. The simplicity of the SPH method allows complex dispersion. The model can be used to simulate complex free-surface flow problems to be solved with relative ease in one, two, flow phenomenon such as invasion of wetting and nonwetting fluids and three dimensions. The SPH was first introduced by into three-dimensional fractures. These processes are a severe challenge for grid-based methods. Surface tension was simulated by using Lucy (1977) and Gingold and Monaghan (1977) for the a van der Waals equation of state and a combination of short-range solution of Navier–Stokes equations in the context of repulsive and longer-range attractive interactions between fluid partiastrophysical fluid dynamics. Since its introduction, SPH cles. The wetting behavior was simulated using similar interactions has been successfully used to model a wide range of between mobile fluid particles and stationary boundary particles. The fluid flow processes and the behavior of solids subjected fracture geometry was generated from self-affine fractal surfaces. The to large deformations. For example, Monaghan (1994) fractal model was based on a large body of experimental work, which used SPH to model the collapse of dams, Morris et al. indicates that fracture surfaces have a self-affine fractal geometry (1997) extended SPH to model low Reynolds number characterized by a material independent (quasi universal) Hurst expoflows, and Zhu et al. (1999) and Zhu and Fox (2001) nent of about 0.75. applied SPH to study pore-scale flow and transport in saturated porous media. Incorporation of the effects of surface tension into SPH simulations has been a vexing U flow in fractures plays an important problem. Morris (2000) modeled surface tension based role in the migration of fluids in the subsurface on its continuous macroscopic description with surface under natural and altered conditions. For example, fractension forces that are proportional to the fluid–fluid tures may significantly decrease the time required for interface curvature. This approach gives an accurate contaminants to migrate through the vadose zone to an estimation of the effects of surface tension but involves underlying aquifer. Computer modeling is playing an rather complex calculations of front curvatures that in ever increasing role in the development of a better unsome cases may lead to significant errors. Nugent and derstanding of fluids and the prediction of fluid flow in Posch (2000) used attractive forces, corresponding to oil reservoirs, hydrothermal reservoirs, contaminated the cohesive pressure in the van der Waals equation of sites, aquifers, and other systems. However, the applicastate to simulate surface tension in two-dimensional SPH tion of standard grid-based numerical methods to freesimulations. They found that it was necessary to increase surface and multiphase fluid flow processes with comthe range of the attractive forces to at least twice the plex dynamical interfaces is fraught with difficulties such range of the SPH weighting function to obtain stable as artificial interface broadening and grid entanglement. vapor drops, a costly solution causing a significant inAn alternative approach is to use particles to represent crease in a computational time. In the work described the fluids. The fluid–air or fluid–fluid interfaces move in this paper we used a different approach. Initially we with the particles, there is no need to explicitly track used van der Waals equation of state to simulate formathe interfaces, and processes such as fluid fragmentation tion and behavior of three-dimensional fluid drops. When and coalescence can be handled without difficulty. we used a low particle density this resulted in the formaIn this paper we describe a numerical model based tion of a liquid drops that did not have a stable smooth spherical shape, a problem similar to that encountered A.M. Tartakovsky, Pacific Northwest National Lab., P.O. Box 999/ by Nugent and Posch (2000). When we increased the MS K1-85, Richland, WA 99352; P. Meakin, Idaho National Lab., particle density (that corresponds to increase in resoluP.O. Box 1625, MS 2025, Idaho Falls, ID 83415-2025. Received 9 Dec. tion of simulation) we were able to create stable spheri2004. *Corresponding author ([email protected]). cal fluid drops but it also led to collapse of some partiPublished in Vadose Zone Journal 4:848–855 (2005). cles. When a combination of short-range repulsive and Technical Reports large-range attractive interactions (with the range of the doi:10.2136/vzj2004.0178 © Soil Science Society of America 677 S. Segoe Rd., Madison, WI 53711 USA Abbreviations: SPH, smoothed particle hydrodynamics 848 Published online August 16, 2005